Transcriptome analysis provides new insights into cold adaptation of corsac fox (Vulpes Corsac)

Abstract Vulpesare widely distributed throughout the world and have undergone drastic physiological and phenotypic changes in response to their environment. However, little is known about the underlying genetic causes of these traits, especially Vulpes corsac. In this study, RNA‐Seq was used to obtain a comprehensive dataset for multiple pooled tissues of corsac fox, and selection analysis of orthologous genes was performed to identify the genes that may be influenced by the low‐temperature environment. More than 6.32 Gb clean reads were obtained and assembled into a total of 173,353 unigenes with an average length of 557 bp for corsac fox. Selective pressure analysis showed that 16 positively selected genes (PSGs) were identified in corsac fox, red fox, and arctic fox. Enrichment analysis of PSGs showed that the LRP11 gene was enriched in several pathways related to the low‐temperature response and might play a key role in response to environmental stimuli of foxes. In addition, several positively selected genes were related to DNA damage repair (ELP2 and CHAF1A), innate immunity (ARRDC4 and S100A12), and the respiratory chain (NDUFA5), and these positively selected genes might play a role in adaptation to harsh wild fox environments. The results of common orthologous gene analysis showed that gene flow or convergent evolution might be an important factor in promoting regional differentiation of foxes. Our study provides a valuable transcriptomic resource for the evolutionary history of the corsac fox and the adaptations to the extreme environments.


| INTRODUC TI ON
Adaptation is a process in which organisms adapt to their environment to survive and reproduce in terms of morphological structures and physiological functions (Brandon, 1978). In particular, creatures in extreme environments have undergone adaptive evolution in phenotype, physiology, and behavior, and the study of its adaptive mechanism has always been a hot topic in field of evolutionary biology. Numerous studies have explored in depth the mechanisms of adaptation in different species at the genomic, metagenomic, and transcriptome levels, such as sheep (Ji et al., 2016), yak (Lan et al., 2021), goat (Martchenko et al., 2020;, and dog (Gou et al., 2014;Miao et al., 2017;Wu et al., 2016).
As a typical example of adaptation, the fox distributed widely and the population remains stable around the world (IUCN, 2021, Figure 1). The members of this genus play important roles in maintaining ecosystems, including tundra, forests, deserts, and grasslands (Statham et al., 2014). Accordingly, they have undergone adaptive evolution in physiology, morphology, and other aspects.
However, little is currently known about their underlying genetic causes, something we aimed to explore in this study. Research on the molecular ecology of the Vulpes has mainly focused on red fox and arctic fox, while corsac fox currently exhibits only in behavior, diet, and zoonoses (Abdybekova et al., 2015;Guilherme et al., 1938;Lkhagvasuren et al., 2016). In the study of adaptive evolution, genome-wide extension of unique gene families and positively selected genes potentially play roles in the adaptation of extreme environments in arctic fox (Peng et al., 2021).
Transcriptome analysis of arctics and the red fox has shown that fatty acid metabolism genes are under positive selection in the arctic fox and seem generally to be crucial for mammalian survival under arctic conditions (Kumar et al., 2015). In our previous study, several positively selected genes related to natural immunity (CFI and LRRFIP1), protein synthesis (GOLGA4, CEP19, and SLC35A2), and DNA damage repair (MDC1) were screened in farmed fox, which might play a role in adaptation to breeding environments (Yang et al., 2021).
Transcriptome analysis is an economic strategy of functional genomic studies that enables the collection of gene sequences for nonmodel animals. These advances greatly facilitate comparative evolutionary studies (Emrich et al., 2007;Garber et al., 2011;Zhang et al., 2010). Despite not representing an entire genome, this transcriptome analysis identified numerous genes that are relevant to adaptation in species. Transcriptome analysis has been used for the assembly and annotation of transcriptomes in many Canidae, including wolves (Liu et al., 2017(Liu et al., , 2019, dogs (Yang et al., 2018), foxes (Kukekova et al., 2011;Kumar et al., 2015;Yang et al., 2021), and Raccoon dogs (Du et al., 2017). These studies have identified many genes involved in adaptive evolution and have contributed to the development of extensive genomic and transcriptomic resources for Canidae. To date, few studies have investigated wild environment adaptation in corsac fox, and the molecular mechanisms have not yet been illuminated. In this study, the pooled transcriptome of multiple tissues of corsac fox was sequenced and compared with red fox and arctic fox to explore the molecular mechanism and corresponding genes of their adaptation to cold environments. The findings provide insights into the molecular mechanisms of adaptation to extreme environments in Vulpes.

| Sample collection
Kidney, brain, and liver tissue of one Vulpes corsac (corsac fox, named CF) was collected from Hulun Lake National Nature Reserve, Inner Mongolia, China. The tissue sample was stored in RNAlater RNA stabilization reagent (QIAGEN, USA) and frozen at −80℃ in ultradeepfreeze equipment. The data of two vulpes lagopus (arctic fox, named AF1 and AF2) and one vulpes vulpes (red fox, named RF) were downloaded from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database with accession numbers ERR687853, ERR687854, and ERR687855, respectively, and this sample was collected from Kronoby, Finland (Kumar et al., 2015).

| RNA extraction, library construction, and sequencing
Total RNA of CF was extracted using RNeasy Mini Kit (QIAGEN, USA) and checked for purity using 1% agarose gels and a NanoPhotometer ® spectrophotometer (IMPLEN, USA). The integrity and concentration of RNA were measured and assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA) and the Qubit ® RNA Assay Kit (Invitrogen, USA) in a Qubit ® 2.0 Flurometer (Life Technologies, USA), respectively. Total RNA (at least 3 mg per tissue) was used to construct sequencing libraries with the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, USA). Once the library was built, an Agilent Bioanalyzer 2100 system was used to assess the F I G U R E 1 Corsac fox in HulunBuir grassland, Inner Mongolia library quality and then sequenced on an Illumina HiSeq 2500 platform.

| Transcriptome assembly, gene functional annotation, and CDS identification
Clean reads were obtained from raw reads by removing reads that contained adapter or poly-N and low-quality reads (Q-value ≤ 20) (Cock et al., 2010). Transcriptome assembly was carried out with Trinity  to obtain reference transcriptome for subsequent analysis. Annotation of gene functional included pathway annotation, protein functions annotation, Gene Ontology (GO) annotation, and COG/KOG functional annotation. The website and parameters of those databases are listed in Table S1. The NR and Swiss-Prot databases were used to blast the CDS of each putative unigene with a cut-off of 1e-5, and Estscan software was used to determine the unigene that did not have aligned results (Iseli et al., 1999).

| Identification of gene orthologous groups and calculation of Ka/Ks
OrthoMCL (Li et al., 2003) software was used to identify orthologous groups with Ailuropoda melanoleuca (aml) and Canis lupus familiaris (cfa) as the outside and internal reference genomes, respectively. Sequence and optimized protein alignment results were constructed by Muscle 3.8.3313 and Gblocks_0.91b (Edgar, 2004 When ω > 1 was usually said to be evolving under positive selection, ω < 1 was usually said to be under purifying selection of homologous genes. In this study, the value of ω was calculated for the pairwise alignments by the Paml-codeml 4.7 (Yang, 2007) package with default settings.

| Homologous gene analysis of three species of fox
The orthologous genes obtained in this study were compared with those screened in previous studies on silver fox and blue fox to acquire the orthologous genes shared by five samples (corsac fox, red fox, arctic fox, blue fox, and silver fox) (Yang et al., 2021). Then, phylogenetic tree was constructed based on orthologous genes and 12 protein-coding genes (except ND6) of mitochondria downloaded from the NCBI database by three methods (PAUP (Swofford, 1998), Bayes (Fredrik & Huelsenbeck, 2003) and RAxML (Alexandros, 2014)), and the reasons for the differences were discussed based on previous studies (Degnan & Rosenberg, 2009;Nakhleh, 2013).

| Overview of transcriptome sequencing data and de novo assembly
One cDNA library was constructed from pooled RNA extracts from different tissues (liver, brain, and kidney) of Vulpes corsac (CF). The summary statistics of the RNA-Seq data are shown in Table 1. After filtering the raw data, a total of 73.37 million clean reads were obtained with an average read length of 150 bp base pairs for further analysis. The Q30% and Q20% of the CF sample were 94.13% and 97.69%, respectively. A total of 234.9 million clean reads remained from RF, AF1, and AF2, which were downloaded from the NCBI database.   Note: N50/N90: The transcript obtained by splicing was arranged from long to short and then accumulated. When the cumulative length>=50%/90% of the total length, then the transcript length is considered N50/N90.

| Calculation of Ka/Ks to identify positively selected genes
To reveal the genetic mechanism of adaptation in response to the environment, Paml-codeml software was used to perform selective stress analysis on the identified single-copy orthologous genes. The 393 identified single-copy orthologous genes were used to calculate the Ka/Ks ratio. The results showed that 16 orthologous genes were identified to be under positive selection (ω > 1, Table S2, Figure 4).

| Homologous gene analysis of three species of fox
Combined with our previously published study, a total of 27 orthologous genes were obtained in five samples. Among them, the  (Table S3, Figure 5b). Figure 5, the mitochondrial phylogenetic tree was the same as the traditional classification, but it was different from the results of orthologous gene construction. That is, the classification status of the corsac fox was different. To address this, gene trees of the 25 orthologous genes were constructed using RAxML and analyzed according to previous studies.

As shown in
A total of 15 gene trees were constructed from 25 orthologous genes ( Figure S1). Among them, the gene trees of ARL3, CCDC151, GPC3, RDH10, and TP53I3 were the same as the tree constructed by mitochondria, and of the rest of the evolutionary tree, there were two main types: (1) the gene tree of 9 genes (ARF4, F13B, OTUD5, RGS4, SCNN1B, ANAPC5, PXN, FAM122A, and CAMLG) was the same as the tree constructed by orthologous genes. That is, corsac fox and arctic fox were the first cluster, with rad fox being the second cluster ( Figure 5c); (2) the gene tree of 6 genes (FAM171B, MGST1, GPLPH3, RARRES2, RND1, and LONRF1) showed that rad fox and arctic fox were the first cluster, with corsac fox being the second cluster ( Figure 5d).

| DISCUSS ION
Vulpes corsac, Vulpes vulpes, and Vulpes lagopus have evolved different physiological, morphological, and behavior adaptations to cope with harsh environments. In the present study, comparative transcriptome analysis was used to identify the selected genes and pathways between them.
The results of orthologous gene enrichment showed that three genes (AQP3, EGR1, and FABP1) were annotated with pathways associated with the hypoxic response. The AQP3 gene encodes aquaporin 3, which is located in the basolateral membrane of renal collecting duct cells and could promote the transport of small non-ionic solutes such as urea and glycerol (Hou et al., 2016).
The FABP1 gene encodes fatty acid-binding proteins, which are highly conserved hydrophobic proteins that are responsible for the uptake, transport, and metabolism of fatty acids (Linssen et al., 1990). In addition, there were also a large number of pathways related to nerve cells, immunity, and protein and lipid metabolism, including a large number of genes. These results suggest that these important homologous genes might have been preserved in response to harsh environments.
Selective pressure analysis results showed that a total of 16 genes were screened under positive selection from 393 single-copy orthologous genes. Enrichment analysis of positive genes showed that the LRP11 gene which encodes low-density lipoprotein receptorassociated protein was enriched in several pathways related to the low-temperature response. In addition to our enrichment results, previous studies have shown that it is involved in animal anxiety responses (Cai, 2015). We speculated that this gene might play a key role in the response to environmental stimuli of animals, and the positive selection of this gene might be the evolution of foxes to adapt to various environmental changes, such as cold and hunger.
In addition, several genes were related to DNA damage repair (ELP2 and CHAF1A), innate immunity (ARRDC4 and S100A12) and the respiratory chain (NDUFA5). Among these genes, the ARRDC4 gene was involved in a pathway involving ubiquitin proteins. This gene has also been shown to promote K63 polyubiquitination of MDA5 through TRIM65 and regulate the innate immune response induced by enterovirus (Meng et al., 2017). The protein encoded by the S100A12 gene acts as an inflammatory mediator molecule that activates RAGE and TLR-4 mediated autoimmune responses. It could also affect the transport of some lipoproteins and increase the calcification of blood vessels (Li & Dewey, 2011). Among the remaining positive selection genes, there were three zinc finger protein genes (ZNF212, KRAB, and ZNF311). The main role of KRAB is to induce gene silencing and maintain genome stability, which plays a huge role in determining the evolution of species (Ivanov et al., 2007;Rowe & Trono, 2011). Studies on orthologous transcription factors in Drosophila have found that there were widespread differences in the binding of this protein to DNA, and such differences occur in a progressive and evolutionarily feasible way (Siggers et al., 2014).
The adaptation of zinc finger proteins combined with the characteristics of new sites enhanced the adaptive evolution of species . These positively selected genes involved in the stress response, DNA damage repair, and innate immunity might play a role in adaptation to harsh wild fox environments.
The results of common orthologous gene analysis showed that the gene tree constructed by orthologous genes was quite different from the species tree. A literature review shows that the main causes of this situation were horizontal gene transfer (HGT), gene duplication and loss, hybridization, and convergent evolution (Degnan & Rosenberg, 2009). Among the possible reasons, HGT is generally believed to occur in microorganisms, such as bacteria, and does not affect phylogenetic relationships in higher animals (Yang, 2017). Sequence alignment results showed that gene duplication or loss did not occur in these orthologous genes. Thus, we speculate that gene flow or convergent evolution might have contributed to this result. (1) If this difference was due to hybridization, then we hypothesize that there was no reproductive isolation between the three species prior to a certain point in time, so gene flow may have existed between them. This gene flow may be an important factor in promoting regional differentiation of species (Yang, 2017).
(2) If these differences were due to convergent evolution, it might mean that these genes evolved convergent in different foxes to adapt to the same environment. In addition, a large number of CDSs discovered in the fox transcriptome data will provide basic transcriptome data for future research.

| CON CLUS IONS
Despite not representing an entire genome, comparative transcriptome analysis was used to identify numerous genes and pathways that are relevant to adaptation in nonmodel organisms. Our study assembled and characterized the first transcriptome data for corsac fox using RNA-Seq technology and revealed candidate LRP11 genes that might be involved in cold environment adaptation in foxes.
The respiratory chain, innate immunity, and DNA damage repair genes under positive selection might play a role in adaptation to

ACK N OWLED G M ENTS
This work was supported by the National Natural Science Foundation of China (Approval No. 31872242, 32070405, 31900311, 32000291, and 32170530).

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflicts of interests.

E TH I C A L A PPROVA L
All sample procedures and experimental methods were approved by the Qufu Normal University Institutional Animal Care and Use Committee (No. QFNU2017-024), Qufu, China.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data analyses during this study are deposited in the Sequence Read Archive database with accession number of SRR16229812.